% Calculate number of applicants of each type by pp level------------------

kappa0 = zeros(Nu,K,tau);

for idx1 = 1:tau
    for idx2 = 1:Nu
        for idx3 = 1:K
            kappa0(idx2,idx3,idx1) = zeta0{idx1}(idx3,1,idx2)*pi(idx1)/...
                (zeta0{1}(idx3,1,idx2)*pi(1) + zeta0{2}(idx3,1,idx2)*pi(2) + zeta0{3}(idx3,1,idx2)*pi(3));
        end
    end
end

for idx1 = 1:tau
    for idx2 = 1:Nu
        for idx3 = 1:K
            Npbytype0(idx2,idx3,idx1) = Np0(idx3)*kappa0(idx2,idx3,idx1)*PTC(idx2);
        end
    end
end

Npbytype0 = sum(Npbytype0)

kappa1 = zeros(Nu,K,tau);

for idx1 = 1:tau
    for idx2 = 1:Nu
        for idx3 = 1:K
            kappa1(idx2,idx3,idx1) = zeta1{idx1}(idx3,1,idx2)*pi(idx1)/...
                (zeta1{1}(idx3,1,idx2)*pi(1) + zeta1{2}(idx3,1,idx2)*pi(2) + zeta1{3}(idx3,1,idx2)*pi(3));
        end
    end
end

for idx1 = 1:tau
    for idx2 = 1:Nu
        for idx3 = 1:K
            Npbytype1(idx2,idx3,idx1) = Np1(idx3)*kappa1(idx2,idx3,idx1)*PTC(idx2);
        end
    end
end

Npbytype1 = sum(Npbytype1)


% Calculate expected opportunity cost by type------------------------------
% S0 = R0_out/abs(mu/200);
% 
% S0_1 = S0(:,2:end,:);
% S0_1(:,10,:) = S0_1(:,9,:);
% 
% for idx1 = 1:tau
%     for idx2 = 1:Nu
% %         OC0(idx2,idx1) = rho*(S0_1(idx2,:,idx1) - S0(idx2,:,idx1))*zeta0{idx1}(:,:,idx2)*PTC(idx2);
%         OC0(idx2,idx1) = rho*(S0_1(idx2,:,idx1) - S0(1,:,idx1))*zeta0{idx1}(:,:,idx2)*PTC(idx2);
%     end
% end
% OC0 = sum(OC0)
% 
% 
% S1 = R1_out/abs(mu/200);
% 
% S1_1 = S1(:,2:end,:);
% S1_1(:,10,:) = S1_1(:,9,:);
% 
% for idx1 = 1:tau
%     for idx2 = 1:Nu
% %         OC1(idx2,idx1) = rho*(S1_1(idx2,:,idx1) - S1(idx2,:,idx1))*zeta1{idx1}(:,:,idx2)*PTC(idx2);
%         OC1(idx2,idx1) = rho*(S1_1(idx2,:,idx1) - S1(1,:,idx1))*zeta1{idx1}(:,:,idx2)*PTC(idx2);
%     end
% end
% OC1 = sum(OC1)

OC0 = zeros(Nu,K,tau);
for idx1 = 1:tau
    for idx2 = 1:Nu
%         OC0(idx2,:,idx1) = sum(P0{idx1}(2:end,:,idx2).*(Vb0{idx1}(2:end,:,idx2) - repmat(Vb0{idx1}(1,:,idx2),[22,1])))/abs(mu/200)*zeta0{idx1}(:,:,idx2)*PTC(idx2); 
        OC0(idx2,:,idx1) = 100*(sum(P0{idx1}(2:end,:,idx2).*(Vb0{idx1}(2:end,:,idx2)./repmat(Vb0{idx1}(1,:,idx2),[22,1])))-1).*kappa0(idx2,:,idx1)*PTC(idx2);
    end
end



OC1 = zeros(Nu,K,tau);
for idx1 = 1:tau
    for idx2 = 1:Nu
%         OC1(idx2,:,idx1) = sum(P1{idx1}(2:end,:,idx2).*(Vb1{idx1}(2:end,:,idx2) - repmat(Vb1{idx1}(1,:,idx2),[22-size(drop,2),1])))/abs(mu/200)*zeta1{idx1}(:,:,idx2)*PTC(idx2);
        OC1(idx2,:,idx1) = 100*(sum(P1{idx1}(2:end,:,idx2).*(Vb1{idx1}(2:end,:,idx2)./repmat(Vb1{idx1}(1,:,idx2),[22-size(drop,2),1])))-1).*kappa1(idx2,:,idx1)*PTC(idx2);
%         OC1(idx2,:,idx1) = (sum(P1{idx1}(2:end,:,idx2).*(Vb1{idx1}(2:end,:,idx2) - repmat(Vb1{idx1}(1,:,idx2),[22-size(drop,2),1])))-1).*kappa1(idx2,:,idx1)*PTC(idx2);
    end
end

sum(sum(OC1)) 

sum(sum(OC0))


% Calculate mean marginal willingness to pay (for Buschena comparison)-----
S = R0_out/abs(mu/200);

S_1 = S(:,2:end,:);
S_1(:,10,:) = S_1(:,9,:);

S0 = repmat(S(:,1,:),[1,K,tau]);

d = 0;

% Calculate denominator of mWTP expression
for idx4 = 1:22
for idx1 = 1:Nu
    for idx2 = 1:tau
        for idx3 = 1:K
            denom{idx4}(idx1,idx3,idx2) = phi0(idx4+1,idx3)*P0{idx2}(idx4+1,idx3,idx1)...
                *zeta0{idx2}(idx3,1,idx1)*pi(idx2)*PTC(idx1);
        end
    end
end
denom{idx4} = sum(sum(sum(denom{idx4})));
end

% Calculate numerator of mWTP expression
for idx4 = 1:22
for idx1 = 1:Nu
    for idx2 = 1:tau
        for idx3 = 1:K
            num{idx4}(idx1,idx3,idx2) = rho*(S_1(idx1,idx3,idx2)-S0(idx1,idx3,idx2))...
                *phi0(idx4+1,idx3)*P0{idx2}(idx4+1,idx3,idx1)*zeta0{idx2}(idx3,1,idx1)*pi(idx2)*PTC(idx1);
        end
    end
end
num{idx4} = sum(sum(sum(num{idx4})));
end

mWTP = zeros(22,1);
for idx1 = 1:22
    mWTP(idx1) = num{idx1}/denom{idx1};
end